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ABSTRACT 

We formulate the problem of the formation and subsequent evolution of fragments (or cores) in magnetically- 
supported, self-gravitating molecular clouds in two spatial dimensions. The six-fluid (neutrals, electrons, 
OO ' molecular and atomic ions, positively-charged, negatively-charged, and neutral grains) physical system is gov- 

, erned by the radiative, nonideal magnetohydrodynamic (RMHD) equations. The magnetic flux is not assumed 

to be frozen in any of the charged species. Its evolution is determined by a newly-derived generalized Ohm's 
law, which accounts for the contributions of both elastic and inelastic collisions to ambipolar diffusion and 
CJ . Ohmic dissipation. The species abundances are calculated using an extensive chemical-equilibrium network. 

2^ ' Both MRN and uniform grain size distributions are considered. The thermal evolution of the protostellar core 

and its effect on the dynamics are followed by employing the grey flux-limited diffusion approximation. Real- 
istic temperature-dependent grain opacities are used that account for a variety of grain compositions. We have 
augmented the publicly-available Zeus-MP code to take into consideration all these effects and have modified 
several of its algorithms to improve convergence, accuracy and efficiency. Results of magnetic star formation 
. simulations that accurately track the evolution of a protostellar fragment from a density ~ 10^ cm^^ to a den- 

Dh' sity ~ lO^^ cm while rigorously accounting for both nonideal MHD processes and radiative transfer, are 

Q \ presented in a separate paper 

' Subject headings: ISM: clouds — magnetic fields — MHD — stars: formation — radiative transfer — dust, 

' extinction 

1. INTRODUCTION — BACKGROUND 

The formulation of a theory of star formation is a formidable task. It requires understanding of the nonlinear interactions among 
self-gravity, magnetic fields, rotation, chemistry (including grain effects), turbulence, and radiation. Stars form in fragments 
within interstellar molecular clouds, which have sizes r anging frorn 1 to 5 pc, masses from a few tens to lO'"^ M0, mean densities 
^ ^ \ ~ 10^ cm""^, and temperatures ~ 10 K (Myers 1985; Heileslll987h . Their spectral lines have Doppler-broadened linewidths that 
suggest supersonic (but subAlfvenic) internal motions. In the deep interiors of such clouds, high-energy cosmic rays (> 100 
MeV) maintain a degree of ionization Xi < 10^^, whereas ultraviolet (UV) ionization is responsible for a much greater degree of 
\ ionization .Tj > 10^^ in the outer envelopes. 
00 . 

' 1.1. Magnetic Fields and Star Formation 

■ The possible importance of magnetic fields to the support of interstellar clouds and to the regulation of star formation was 
l>( first studied by Chandrasekhar & Fermi (1953), Mestel & Spitzer (1956), and Mestel (1965) using the virial theorem. Simi- 

lar investigations by iStrittmatteil (I1966aibl) and iSpitzerl (119681) followed. iMestell (|1966l|) calcul ated the magnetic forces on a 

■ spherically-symmetric, gravitationally-bound cloud. Self-consistent calculations by lMouschoviasI (Il976a lbl) produced exact equi- 
■ ■ ' Ubria of initially uniform, isothermal, magnetic clouds embedded in a hot and tenuous, electrically-conducting external medium. 

iMouschovias & Spitzerl (Il976h used these equilibrium states to find the critical mass-to-flux ratio (Af/<I'B)cr = 1/ (BSG)^/^ that 
must be exceeded for collapse against the magnetic forces to set in. Scott & Black (1980) performed numerical simulations of 
the collapse of a supercritical (as a whole) m agnetic cloud. A pic ture of molecular clouds emerged in which magnetic fields play 
a central role in their support and evolution dMouschoviasll 1 978i) . To this date, magnetic braking remains the only mechanism 
that has been shown quantitatively to resolve the most significant dynamical problem of star formation, namely, the angular 
momentum problem (Mouschovias & Paleologou 1979, 1980; see also summary below). 

Subsequent observations lent credence to this picture by revealing the importance of magnetic fields through both dust polariza- 
tion measurements and Zeeman observations. Pol arization studies have exhibited large-scale ordered magnetic fields connecting 
proto stellar cores to their surrounding envelopes dVrba etalJI 19811; iHever et alJll987H Novak et al."1989l, '1997'; 'Lai et al ' 
2003t IC rutcher et al. 2004; Alves et al. 20081 ), often with an hourglass morphology jSc hleuning 1998; Hildebrand et al. 
Girart e t al. 1999; Schleuning et al. 2000; Lai et al ."2002'; 'Cortes & Crutcher'2006'; Girart et al.ii2006) . as predicted by the theoret- 
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ical calculations (Mouschovias 1976b; Fielder & Mouschovias 1993). A large body of Zeeman observations (Crutcher & Kazesj 
1^83; Kazes & Crutcher 1986; Troland et al. 1986; Crutcher et al. 1987; Goodman et al. 1989; Crutcher et al. 19931 1199411 9961 
Troland et al... 1996:, Crutcher et aL,1999a.b : Crutcher 1999: Heiles & Crutcher..2005; Cortes et al...2005.) revealed magnetic fields 



in the range ~ 10 — 200 /j,G in molecular clouds, from small isolated ones to massive star-forming ones. These values are more 
than sufficient to establish the importance of magnetic fields in molecular cloud dynamics. 

It was recognized early on (e.g., see Babcock & Cowling 1953, p. 373) that the magnetic flux of an interstellar blob of mass 
comparable to a stellar mass is typically several orders of magnitude greater than that of magnetic young stars. This is the so-called 
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"magnetic flux problem" of star formation. It lies in the fact that substantial flux loss must take place at some stage during star 
formation. Ambipolar diffusion (the relative motion between plasma and neutrals) was first proposed by Mestel & Scitzer ( 195^ 
as a means by which an interstellar cloud as a whole would reduce its magnetic flux and thereby col lapse. iPneuman & Mitchelll 
(|l965) undertook a detailed calculation of the collapse of such (spherical) cloud. ISpitzerl (Il968h calculated the ambipolar- 
diffusion timescale by assuming that the magnetic force on the ions is balanced by the (self-)gravitational force on the neutrals. 
iNakanol (fT979) followed the quasistatic contraction of a cloud due to ambipolar diffusion using a sequence of Mouschovias' 
(1976b) equilibrium states, each one of which had a smaller mag netic flux than the previous one. 

A new solution for ambipolar diffusion b v iMouschoviasI (1 1 979i) showed that the essence of this process is a redistribution of 
mass in the central flux tubes of a molecular cloud, rather than a loss of magnetic flux by the cloud as a whole. He found the 
ambipolar-diffusion timescale to be typically three orders of magnitude smaller in the interior of a cloud than in the outermost 
envelope, where there is a much better coupling between neutral particles and the magnetic field because of the much greater 
degree of ionization. This suggested naturally a self-initiated fragmentation of (or core formation in) molecular clouds on the 
ambipolar-diffusion timescale 

TAD = 1.8 X 10^ (^) yr 

(Mouschovias 1979, eq. [37]). The inefficiency of star formation was thereby attributed to the self-initiated formation and 
contraction of molecular cloud fragme nts (or cores) due to ambipolar diffusion in otherwise magnetically supported clouds 
(IMouschoviasI 1 976R ll977Lll978l fl979l) . The central mass-to-flux ratio eventually exceeds its critical value for collapse. 
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(see Mouschovias 1976a, eq. [44]), and dynamic contraction ensues. Detailed numerical calculations in slab 
dPaleologou & Mouschovias| |T983; Mouschovias et al. 1985 J), cyUndrical (Mouschovias & Morton 1991, 1992a. b|), and axisym- 
metric geometries (Fielder & Mouschovias 1992, 1993; Ciolek & Mouschovias 1993, 1994, 1995; Basu & Mouschoviaslll994 
Il995afb ; Ciolek & Konigl 1998; Desch & Mouschovias 2001; Tassis & Mouschovias 2005ab, 2007a,b.c) transformed this sce- 
nario of star formation into a theory with predictive power. 

1.2. Rotation 

During the early, isothermal phase of star formation, a cloud (or a fragment) must also lose a large fraction of its angular 
momentum (e.g., see Spitzer 1968, p. 231). Observations show that molec ular clouds and embedded fra gments (or cores) rarely 
exhibit rotation significantly greater than that of the background medium (iGoldsmith & Arquillalll985h . If angular momentum 
were conserved from the initial galactic rotation (i.e., starting from angular velocity fio — 10~^^ s~^ at the mean density of 
the interstellar medium ~ 1 cm"'^), centrifugal forces would not allow even the formation of interstellar clouds (Mouschovias 
1991, § 2). Fragmentation does not alter this conclusion (Mouschovias 1977, § 1). This is referred to as the "angular momentum 
problem" of star formation. As far as clouds and their cores are concerned, the angular momentum problem has been shown to 
be resolved by magnetic braking (i.e., the transport of angular momentum from a fragment to its surrounding medium through 
th e propagation of torsional Alfven waves a long magnetic field lines connecting the fragment to the cloud e nvelope) analytically 
bv lMouschovias & Paleologoul ( II979llI980l) and numericallv bv Bas u & Mouschov ias ( 1994. 1995a b) and lMellon & Lil ( |2008|) . 
Hence the centrifugal forces resulting from the cloud's or core's rotation have a negligible effect on the evolution of the contracting 
core, at least up to central densities of « 10^''^ cm^'^ (see the last paragraph of Tassis & Mouschovias 2007b). 

1.3. Grain Effects 

Inte rstellar grains comprise about 1% of the mass in the interstellar medium (IS pitzeijri978h . |Bakeij ( |1979t ) and lElmegreenI 
(|1979) sug gested that charge d gr ains may couple to the magnet ic field and thereby play a role in ambipolar diffusion and star 
formation. lElmegreenI (1 19791) and lNakano & Umebayashil (Il980() compared and ambipolar-diffusion timescale and the free-fall 
timescale and concluded that ambipolar diffusion occurs over too long a tim escale (roughly 10 times greater than free-fall) 
to be a relevant process in star formation. Refinements by the same authors (Elmegreenl l 19861 : lUmebayashi & NakanollT990l : 
iNishi et al. 1991) led to similar conclusions. Thr ough detailed num erical simulations of core formation and evolution including 
the effects of (negative and neutral) dust grains, ICiolek & MouschoviasI (fl993l 11994 found that grains lengthen the timescale 
for the formation of a core because of grain-neutral collisions, but cautioned that the ambipolar-diffusion timescale should not 
be compared to the free -fall timescale in determining its relevance in magnetically-supported clouds, as originally pointed out 
by iMou schovias' (1977), because molecular clouds are not free-falling. Velocities characteristic of such collapse have not been 
observed. Ciolek & Mouschovias ( 1995) extended these calculations by including UV ionization and a variety of atomic metal 
ions (C+, S+, Si+, Mg+, Na+, F e"*"). Attention was also paid to the complementary effect of protostellar evolution on the 
microscopic physics and chemistry (ICiolek & Mouschoviasll996lll998l) . 

1.4. MHD Waves and/or Turbulence 

The extent to which turbulence (or waves) may or may not affect the e volution of a protostellar fragment has been a t opic of 
debate for several decades, receiving increased attention in recent years (iMac Low & Klessenll2004l : IMouschovias et alJ l2006). 
A consensus has yet to be reached concerning even what causes and maintains the observed broad linewidths long thought to be 
indicative of supersonic turbulence ( Zuckerman & Evans 1974; Arons & Max 1975; Larson 1981; Zweibel & Josafatsson 1983]; 
iMouschoviai 19871: lMouschovias & Psaltislll995l: IMvers & Gammidl 19991: lElmegreen & Scaloll2004): IMouschovias et al.ll200d) . 
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althoughl Mouschovias & PsaltisI d 19951) and lMouschovias et alj ( l2006l) showed quantitatively that observations contradict the key 
assumption of turbulence simulations, namely, that molecula r clouds are magne tically supercritical by a factor 4 — 10. Despite a 
lack of agreement on the origin of the linewidths, analytical ('M ouschoviasll 1 99 li) and numerical calculations (Eng 2002; Eng & 
Mouschovias 2009, in preparation) have demonstrated that turbulence (or waves) plays an insignificant role in the star formation 
process once dynamical contraction of a fragment (or core) ensues. Observations showing nar rowing and eventual thermalization 
of linewidths in protostellar cores (iBaudry et al.lll98ll Mye rs & Benson 1983[: [ Myers et al.lil983; .Bacmann et al.i,2000.) are in 
agreement with this conclusion and with an earher version of it dMouschoviasfl 9871) ! 

1.5. Radiative Transfer 

During the early phases of star formation the energy produced by compressional heating is radiated away by the dust grains in 
the infrared. At higher densities (~ 10^^ cm^"^), the core traps and retains part of this heat and its temperature begins to rise. The 
evolution of the temperatu re in this nonisothermal regime may be approximated (but substantiall y overestimated) by using an 
adiabatic equation of state (lBosslll98lllDesch & M ouschovias"2 001l:lTassis & Mouschoviasll2007al[bll3) . More realistic equations 
of state have also been employed by, for example. Bate (1998). To accurately model the nonisothermal phase of protostellar 
contraction, however, one needs to include a proper treatment of radiative transfer Early efforts to include radiative tran sfer in 
(nonmagnetic) star formation calculations were confined to the use of the diffusion appro ximation (Bodenheimer 196 8; LarsonI 
[l969l[T97llBlack & BodenheimedlT97l IT976t lTscharnute3IT97l lYorke & Krugell[T977h . While the diffusion approximation is 
strictly applicable only to optically thick regions, its ease of implementation and relatively low computational cost make it an 
attractive choice. The Eddington approximation offers a slight improvement in that it retains some of the rigor of using moments 
of the radiative transfer equation, while making the simplifying assumption that the radiation field is ever ywhere isotropic. Its use 
in numerical calculations of (nonma gn etic) star formation has been documented in.Tscharnuter ( 1978|) , iTscharnuter & Winkler! 
(ll979l) . IWinkler & Newmanl(ll980allM . lBossl(ll984ll986ill988l) . and lBoss & Mvhilli(ll995l) . Bvimplicitlv assuming that photons 
always travel a distance comparable to their mean-free path (even if this distance exceeds the free-flight distance cAt, where At 
is the computational timestep), the Eddington approximation gives unphysical behavior in optically thin regions, in which the 
mean-free-path is huge. The result is a signal speed unbound by the speed of light, i.e., it violates causality (e.g., see § 97 of 
Mihalas & Mihalas 1984). 

Increasing the accuracy and realism of a radiative transfer algorithm often requires making limiting assumptions about the 
hydrodynamics in order to make the problem tractable (e.g., Yorke 1980; Masunaga et al. 1998). A full frequency- and 
angle-dependent treatment of the radiation is nearly always confined to postprocessing the results of a hydrodynamic calculation 
(l%rkdll977t lYorke & Shustov 1981; Adams & Shu 1985, 1986) or a grey (i.e., independent of frequency) radiation hydrody- 
namic calculation (Boss & Yorke 1990; Bodenheimer et al. 1990). By contrast, the flux-limited diffusion (FLD) approximation 
dLevermore & Pomra ning 1981) is a propitious compromise that retains some of the advantages of the diffusion and Eddington 
approximations, while preserving causality and coupling self-consistently to the hydrodynamic equations. 

1.6. Outline 

In this paper we formulate the problem of the formation and evolution of protostellar fragments (or cores) in magnetically- 
supported, self-gravitating molecular clouds, including the effects of both ambipolar diffusion and Ohmic dissipation (which 
becomes important at high densities), grain che mistry and dynam ics, and radiation. Using the results of Eng (2002) and Eng 
& Mouschovias (2009, in preparation), and Bas u & MouschoviasI ([l994), we may safely ignore the effects of turbulence and 
rotation, respectively, on the evolution of the protostellar core for the densities considered here. The physical and chemical 
properties of the model cloud are summarized in Section |2] The radiation magnetohydrodynamic (RMHD) equations governing 
the evolution of the model cloud are presented and discussed in Section |3] In Section |4] we present the chemical model used in 
the calculations. The physics of magnetic diffusion (ambipolar and ohmic) is handled by using a generalized Ohm's law, which is 
derived in Section |5] We treat the radiative transfer using the grey FLD approximation, with realistic grain opacities accounting 
for a variety of grain compositions (§ |6]l. The numerical method of solution is discussed in Section |7] Finally, we give the 
simplified set of equations and a brief summary in Section|8] Details, mostly mathematical, are left for the Appendix. Results are 
presented in a separate paper (Kunz & Mouschovias 2009, in preparation). 

2. BASIC PROPERTIES OF THE MODEL CLOUD 

We consider a self-gravitating, magnetic, weakly-ionized model molecular cloud consisting of neutral particles (H2 with 20% 
He by number), ions (both molecular HCO+ and atomic Na+, Mg+, K+), electrons, singly negatively-charged grains, singly 
positively-charged grains, and neutral grains. Following iDesch & Mousc hovias (2001), the abundances of all species (except 
the neutrals) are determined from the chemical reaction network shown in Table [T] and described below in Section |43] Cosmic 
rays of energy > 100 MeV are mainly responsible for the degree of ionization in the cloud. Once column densities ^ 100 g 
cm~^ are achieved, cosmic rays are appreciably attenuated. At even higher densities, cosmic rays are effectively shielded and 
radioactive decays become the dominant source of ionization. Finally, at temperatures on the order of 1000 K or higher, thermal 
ionization of potassium becomes impor tant. UV radiation provides an additional io nization mechanism, but it only affects the 
outer envelope of molecular clouds (HoUenbach et al.ll971tlGrassgold & Langedl974} . We consider spherical grains whose radii 
are determined by either a uniform or an MRN dMathis etal.lll977h size distribution. In the case of collisions of ions (molecular 
or atomic) with grains, we assume that the ions do not get attached to the grains, but rather that they get neutralized, with the 
resulting neutral particle escaping into the gas phase. Thus the total abundance of metals as well as the total HCO abundance 
remain constant. Grain growth is not considered here. 

We follow the ambipolar-diffusion-initiated evolution of an axisymmetric two-dimensional model cloud from typical mean 
molecular cloud densities (~ 300 cm^'^) to densities characteristic of the formation of a hydrostatic protostellar core. The axis of 



4 



ti2 + OK 


H„ +e 

1U"T" 1 T T 

Hg +H 


1 U 

H-o + n2 
H+ + CO 


-+ HC0+ + H2 


HCO+ + e - 


H + CO 


A+ +e 


A" + hu 


A0 + HCO+ - 


^ A+ + HCO 


e + go 


g- 


e + g+ 


go 


A+ + g- 


AO+go 


A+ + go 


AO + g+ 


HCO+ + go - 


HCO + g+ 


HCO+ + g_ - 


HCO + go 


S% + S-' 




- So+go' 


g± + go ' 


go + g± 



TABLE 1 

Chemical reaction network used in the calculation of the abundances of charged species. 

Relevant Chemical Reactions in Molecular Clouds 
Cosmic-Ray Ionization: 

Dissociative Recombination: 
Radiative recombination: " 
Charge transfer: " 
e~ attachment onto grains: 

Atomic-ion attachment onto grains: " 

Molecular-ion attachment onto grains: 

Charge transfer by grain-grain collisions: 

" A"'" represents an atomic ion, such as Na"'" , Mg+ , and K+ , and A" the corresponding neutral atom. 

symmetry is aligned with the z-axis of a cyUndrical polar coordinate system (r, (j>, z). Isothermality is an exc ellent app roximatio n 
for the early stages of star formation, while the density is smaller than w 10 cm -3 (Gaustad 1963-|Hayashi! 1966; La rsonI 19691) . 
However, once the heat generated by released gravitational energy during core collapse is unable to escape freely (at a central 
number density of nopq ~ 10^ cm~^), radiative transfer calculations are employed to determine the thermal evolution of the core. ' 
This is an improvement over previous magnetic star formation calculations to reach these densities (iDesch & Mouschoviasl200"l1: 
ITassis & Mouschovias 2007a, b.c), which assumed an adiabatic equation of state beyond a critical density because of the high 
computational expense of radiative transfer calculations. Numerical techniques and computer hardware have matured enough by 
now to render these once impractical calculations feasible. 

3. THE SIX-FLUID RMHD DESCRIPTION OF MAGNETIC STAR FORMATION 

The RMHD equations governing the behavior of the six-fluid system (neutrals, electrons, ions, negative, positive, and neutral 
grains) are: 

% + V.(p„t;„) = 0, (la) 

+ V . {p,_v,_ + p,,v,„ + p,,v,,) = , (lb) 

^%7^ + V • [pM^) = - VPn - pVi^ + -jxB + -xjrT , (Ic) 
ot c c 

= -en^(^E+^xB^ +F^n, (Id) 
= +eni (E+ — xB]+Fin, (le) 



= -eng_ + ^ X b) + Fg_„ + Fg_g„,inei , (If) 

- +e7ig^ (^E+'^xB^+ Fg^„ + Fg^g„,i„ei , (Ig) 

— -fgon + J^gQg_.incl + gog+,incl j (Ih) 

VxB = —j, (li) 

c 

j = e (fiiVi - TicVc + %+'yg+ - ng_Vg_) , (Ij) 
dB 

— = -cVxE, (Ik) 

= 4^Gp„ , (11) 

+ V • (Un^^n) = -fn V • - 47rKp6 + CK^£ + Tdiff , (Im) 

' We have varied the density riopq at which we turn on the radiative transfer solver from 10® — 10^^ cm"'' and found that Wopq ^ lO'^ cm"'' is necessary 
to achieve a smooth transition from isotheiTnality. This numerical necessity does not mean that the isothermahty assumption breaks down at as low a density as 
10^ cm-3. 
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dE 

— + V • [Ev^) = - V • :F - Vt;„ : P + At^k^B - ckeS , (In) 
ot 

— +V-{J^v,,) = -c^V-P-cxrJ^. (lo) 

The quantities p^, Ug, and Vg refer to the mass density, number density, and velocity of species s; the subscripts n, i, e, g_, g_|_, 
and go refer, respectively, to the neutrals, ions, electrons, negatively-charged grains, positively-charged grains, and neutral grains. 
The quantities E and B denote the electric and magnetic field, respectively, j the total electric current density, Un the internal 
energy density, the gas pressure, and tp the gravitational potential. The source term Fdifi in the internal energy equation (llml i 
represents heating due to ambipolar diffusion and Ohmic dissipation (see § l5.5l below). The magnetic field satisfies the condition 
V • B = everywhere at all times. The definitions and derivations of Fg^ (the frictional force per unit volume on species s 
due to elastic collisions with neutrals) and F^s,inc\ (the force per unit volume o n grain fluid 7 due to the conversion of dust 
particles of fluid S into dust p artic les of fluid 7) are di scussed in detail in § 3 of iTassis & MouschoviasI (l2005a) . as well as in 
ICiolek & MouschoviasI d 19931) and iMouschoviasI (Il996h (in the absence of positively-charged grains). 

The radiation variables are the Planck function B, the total (frequency-integrated) radiation energy density £, the total 
(frequency-integrated) radiation momentum density JF, and the total (frequency-integrated) radiation pressure tensor P: 

£{x,t)^-J dv (j) dni{x,t;n,u) , (2a) 

{x,t) = J dv j> dni{x,t]n,v)n, (2b) 
I f 

P(x,t) ^ — / dv i) dVl Iix,t;n,v) hfi . (2c) 
c Jo / 

Here we have introduced the frequency v, the extinction coefficient (i.e., opacity) x{^) [= + f (i^), where n is the absorption 
coefficient and a the scattering coefficient], and the radiation specific intensity /. The material properties Kp, ke, and xi^ are the 
Planck and energy mean absorption coefficients, and the flux-weighted mean opacity, respectively; they are given by 

Kp = —J K{v)B{v)dv , KE=-J K,{v)£{v)dv , xr=-^j x{v)^{v)dv . (3a,b,c) 

Equations (fTnl i and ( [Tol l are obtained from taking moments of the radiation transport equation 

Id \ 

-— + • V j I{x, t; fl, v) = x{x, t; fl, v) [S{x, t; VI, v) - I{x, t; VI, u)] (4) 

under the assumptions that all the radiation variables are measured in the comoving frame of the fluid (in this frame the material 
properties are isotropic) and that the material properties are grey (Mihalas & Mi halas 1984,) . ^ We have taken the source function 
S in the transport equation (|4]i to be given by 

AirS^ = ■ , (5) 

taking into account both establishment of local thermodynamic equilibrium and coherent isotropic scattering of radiation 
(iMihalas & MihalaslfT98l . 

In the force equations for the electrons, ions, and grains, the accel eration terms have been ne glected due to the small inertia of 
these species. The acceleration term for the plasma was included bv' Mouschovias et al.l (11985 ') and it was shown that the plasma 
reaches a terminal drift velocity very fast. Similarly, the thermal-pressure and gravitational forces have been dropped from the 
force equations of all species other than the neutrals because they are negligible compared to the electromagnetic and collisional 
forces. The inelastic momentum transfer by the electron and ion fluids due to attachment onto grains and neutralization are 
negligible compared to the momentum transfer due to elastic collisions, and they have been omitted from the force equations ( fTdl i 
and ( [Tel l (see discussion in § 3.1 of Ciolek & Mouschovias 1993). As we consider a distribution of grain sizes, it should be noted 
that equations ( [Tbb . (fTfl) - (fThl l apply to each grain size separately. 

The full set of RMHD equations are closed with constitutive relations for the gas pressure, opacities, and the Planck function 
[i.e., Pn ~ Pn{pn,T), XT = XJ^iPg^T), ke = KE{pg,T), Kp = Kp{pg,T), and B ~ B{T), where T is the gas temperature 
and pg = pg_ + pg„ + pgj_ is the total grain mass density]. In addition, we close the radiation moment equations with the tensor 
variable Eddington factor f which is used to eliminate the radiation stress tensor P in favor of the radiation energy density £ via 

P = f£ . (6) 



The Eddington factor f is determined by employing the flux-limited diffusion approximation (see § l6.1l below). The equation of 
state for an ideal gas is given by 

P„ = (7 - , (7) 



^ We caution here that Preibisch et al] 09951) and|Yorke & Sonnhalter (2002) have shown that multi-frequency calculations generally produce higher dust 
temperatures and greater degrees of anisotropy in the radiation field than corresponding grey calculations. 
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Fig. 1 . — Dependence of 7 (ratio of specific heats) on temperature for H2. 

where 7 = fee /cv + 1 is the adiabatic index and cy is the specific heat at constant volume per particle: 

CM^\k^+ct + d^\ (8) 

assuming that there is no coupling between the rotational and vibrational degrees of freedom of the molecule (in this case, H2). 
The vibrational specific heat is 



SvibV exp(evib/7') 



(9) 



Cv* = TfcBX^^lnZo + -fcBx2£-lnZp, (10) 



T J [exp(evib/T)-l]' 
where O^ib = 6100 K; the rotational specific heat for a 3:1 mixture of ortho- and para-hydrogen is 

3 , 2 5' , ^ 1 , 2 5' 

where Zq and Zp are the ortho- and para-hydrogen partition functions, respectively, given by 

Zo=Yl (2i + 1) exphxj(j + 1)] and = (^J + 1) exp[-a;j(j + 1)] , (1 la,b) 

odd J even J 

and X = Q^ot/T — 85.4 K/T (lKittellll938h . The dependence of 7 on temperature T is shown in Figure[T] 

Altogether, then, we have a system of 17 equations [( [Tal l - ( fTob . and (|7])], which contain 21 unknowns {pn, Pn, u^, E, B, 
j, tp, Vn, Ve, Vi, Vg_, Vg^, Vgg, pc, Pi, Pg_, Pgj_, Pgo , E, P). To closc the system, the densities of electrons, ions, and charged 
grains (n^, ni, ng_ , and ng_^) are calculated from the equilibrium chemical model described below. 

4. THE CHEMICAL MODEL 

4.L Ionization Rate 

The rate of ionization per unit volume is given by C,nYi_^ and is (in principle) due to the following ionization sources: ultraviolet 
radiation, cosmic rays, radioactivities, and thermal ionization. The presen ce of molecules in molecular c louds implies low levels 
of UV radiation, so it is usually neglected. UV radiation was included bv lCiolek & MouschoviasI ( Il995h in their numerical sim- 
ulations of core formation and evolution. They found that UV ionization dominates cosmic-ray ionization for visual extinctions 
Ay < 10 and can increase the degree of ionization in the envelope by at least two orders of magnitude (see also McKee 1989). 
The increase in ionization was found to speed up core collapse by approximately 30% because the central gravitational field of a 
flattened cloud is stronger (i.e., less diluted by the mass of the envelope) when matter in the envelope is held farther away from 
the forming core. Once dynamical contraction ensues, it was found that UV radiation has little effect on the evolution of central 
quantities and therefore it is usually neglected. For numerical reasons, however, we add to the electron and ion number densities 
the second term in equation (4h) of iFielder & Mousc hoviaH (Il992|) (— 467.64 njj^ cm^'^) so as to maintain a relatively large 

degree of ionization (^ 10^^ — 10~^) (and therefore negligible ambipolar diffusion) in the low-density (nng S 10"^ cm^"^) cloud 
envelope. This term qualitatively mimics the effect of cloud envelope penetration by UV photons, and has negligible quantitative 
effect on the formation and evolution of the core. 

Cosmic rays, on the other hand, with typical energies of 100 MeV are able to penetrate deeper into molecular clouds. 
lUmebavashi & Nakanol (Il980h have investigated the ionization due to a spectrum of cosmic rays at various energies. They 
found that the cosmic-ray ionization rate was well described by the following relation: 

CcR = Coexp(-SHj96gcm-2), (12) 
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(nK+) = 4.1 X 10 l^TlH^^K" innnr^ Cm-* « . (13) 



where is the column density of H2 separating the point in ques ti on from the exterior of t he clou d and (^0 = 5 x 10^^^ s^^ 
is the canonical unshielded cosmic-ray ionization rate (ISpitzerlll978l) . iTassis & MouschoviasI ( l2007bl) found that, when a typical 
core's central density exceeds ~ 10^^ cm^"^, cosmic rays are shielded and an abrupt decrease in ionization occurs. 

Once the core is shielded from high-energy cosmic rays, the dominant source of ionization is radioactive decay of "'"K or 
^® Al. The isotope ^"K is the most common radionuclide invoked, due to its long half-life of 1 .25 Gyr and its ubiquity in nature 
(0.012% of terrestrial potassium is ''°K). The density of potassium in the interstellar medium (2.70 x 10^^ n-Ha) and the energy 
of the beta particle emitte d as "'^°K decays, 1.31 MeV, are u sed as inputs to calculate the ionization rate (e.g., see Glassgold 1995): 
^40 = 2.43 X 10^^^ s^^. IConsolmagno & Jokipiil (Il978h have suggested that ^^Al may have been a much more potent ionizer 
than ''"K. Performing a similar calculation for ^^Al, one finds (^26 — 1. 94 x 10~^^ s~^, \yith th e fraction of aluminum in the 
isotope ^^Al inferred to have existed in the solar nebula being 5 x 10~^ (Clavton & Leising 1987). Although ^*^A1 is four orders 
of magnitude more potent an ionizer than ^"K, its short half-life (0.716 Myr) makes it relevant only if the initial mass-to-flux 
ratio of the parent cloud is close to critical, so that the evolution is rapid enough to retain an adequate amount of this radionuclide. 
^^Al can also become important if the core happens to get enriched because of a nearby Supernova explosion. 

Finally, at temperatures on the order of 1000 K or higher, collisions between molecules are energetic enough to ionize those 
atoms with low ionization potentials, of which potassium and sodium are the most abundant. The abundance of sodium in the 
interstellar medium is greater than that of potassium (by a factor ~ 14; Lequeux 1975), but the lower threshold of pota ssium 
(4.34 eV vs. 5.13 eV for sodium) makes it the dominant ion. The ionization occurs at a rate jPneuman & Mitchellll 19651) given 
by 

Because this process relies on collisions between two species, it is not expressed in terms of a quantity 

4.2. Grain Size Distribution 

Since the dust opacity, the conductivity of the gas, and the collision rates (see below) all depend on the (local) grain surface 
area, it is necessary to investigate the effect of a grain size distributio n. The initial size d istributions adopted here are a uniform 
distribution and the standard "MRN" distribution of interstellar dust dMathis et alj|1977h . In both cases, the density of the solid 
material of each grain is taken to be pg = 2.3 g cm^'^, the average density of silicates. For the uniform distribution, a fiducial 
grain size gq — 0.0375 fim is used and the total mass density of dust Pg.tot = 0.01pn,tot- For the MRN distribution, the number 
density of spherical dust grains with radii between a and a + rfa is 

dn-g,tot(a) = NuRNa^^-^da . (14) 

The distribution is truncated at a lower grain radius amin and an upper grain radius flmax- The coefficient ^Vmrn is proportional 
to the dust-to-gas mass ratio in the cloud. Note that most of the grain surface area is contributed by small grains, because of their 
overwhelming abundance. 

The grains are binned according to size and charge and treated as separate grain species. Each size bin represents a subset of 
the original distribution of grains, those with radii between aiower and fluppd . The subset of grains in the ath (a = 1,2,..., N) 
size bin is replaced by a number density rtg of grains with identical radii, a^. The total number of grains and the total surface 
area of grains in the size bin are constrained to match the total number and surface area of original grains incorporated into the 
size bin. Hence 

-da, n"a^ = / a — y — da. (15a,b) 



da ^ " -/aiowoi da 

Applying these relations to the MRN grain size distribution, equation ( fT4] l. if there are N size bins, then the ath bin is character- 
ized by grains of number density and radii as follows: 



1 f2.5/Af\ 



f2.5{a-i)/N , 2 ^ 1 a -a ■ r 



0.5/Ar 



1/2 



(16a,b) 



The ratio of the lower and upper radii of the distribution is denoted by ^ = aniin/a,„ax- The total number density of dust, ng^tot, 
is determined by constraining the total grain mass density in the size distribution to be Pg tot: 





[I 


(\- 













eO.5 



(17) 



The lower and upper cutoffs to the size distribution are chosen to be Omin — 0.0181 /im and Omax = 0.9049 [im, respectively. In 
equation ([TtI i. the total mass density of dust in the system, Pg^tot, is chosen in such a way that the total grain surface area in the 
size distribution is equal to that in the fiducial case of a single grain size oq. This constraint demands that Pg.tot be increased by a 
factor (amin/flo) ^^" '"^ over the fiducial value of Pg,tot- Only in this way can the effect of a size distribution, as distinct from just 
the surface area of grains, be determined. With the fiducial values ao = 0.0375 /im and /9g,tot = 0.01/5n,tot for the single grain 
case, Pg.tot — 0.0341pn tot for the case of a size distribution. Empirically, it was found that a minimum of five size bins of grain 
radii were required for convergence of 1%. Since each size grain can be found in one of three possible charge states (— e, 0, and 
+e), a total of fifteen grain species are considered. 

While we do not consider grain growth (and therefore fix the number of grains within each size bin), we do expect the grain 
size distribution to evolve spatially within the star-forming cloud. Ambipolar diffusion can alter a grain size distribution by acting 
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more effectively on the larger grains, causing a spatial segregation of grain sizes that le aves the smaller grains behind i n the cloud 
envelope. The result is a deficit of small grains (a < 10~^ cm) in the cloud core. In fact. lCiolek & Mouschoviasl(ll996l) show how 
observations of grain abundances in the core and envelope of a molecular cloud can, at least in principle, be used to determine 
the initial mass-to-flux ratio of the cloud. 

4.3. Chemical Network 

We use a chemical equilibrium network accounting for electrons (subscript e); molecular ions such as HCO+ (subscript m+); 
neutral metal atoms (subscript A") and atomic ions (subscript A+) of Mg, Na, and K; singly positively-charged grains (subscript 
g+); singly negatively-charged grains (subscript g_); and, neutral grains (subscript go). Multiply negatively- (positively-)charged 
grains may be neglected, because a singly ne gatively- (pos itively-)charged grain repels electrons (ions) thereby decreasing the 
rate of capture by the factor exp(— e^/ ak^T) ( ISpitzeill941l) . The equilibrium assumption is accurate provided that the dynamical 
timescales of interest are sufficiently longer than the chemical-reaction timescales. This is always the case for the density regime 
considered here. The relevant reactions are given below and explained briefly. 

The production of molecular ions (such as HCO+) is balanced by their destruction through charge-exchange reactions with 
atomic ions, by dissociative recombinations (collisions with electrons), or by collisions with and neutralization on the surfaces of 
grains: 

CnH2 = ?^,n+n-At)/3 + nm+neadr + ^?^m+"g°am+g^ + X! "™+"'S?'^ni+g° • (18) 

The index a denotes a grain size bin, and the sum is over all the size bins, which are treated as independent grain species. The 
production of atomic ions by charge-exchange reactions is balanced by radiative recombinations and by collisions with grains: 

n„i+7lAo/3 = TT-A+noarr + ^"A+ng^aA+g" + X^^A+^g^ttA+g" ■ (19) 

a a 

If the atomic ion in question is K+, there is the additional source term due to thermal ionization of potassium atoms, 
nKonH2Q!K"H2- Positively-charged grains are formed by the collisions of ions and neutral grains and by charge exchange be- 
tween grains; they are destroyed by collisions with electrons, collisions with negative grains, and by charge exchange with 
neutral grains: 

m-i- go m-rgg AT go ATgo gc. gf, gQ gQ c g^ cg^ ^ g^ go g°gZ 8+ gg g'^gg ^ ' 

a' a' a.' 

Here the index a' runs over all the grain size bins, independently of the index a. Negatively-charged grains are formed by the 
collisions of electrons and neutral grains and by charge exchange between grains, and are destroyed by collisions with ions, 
collisions with positive grains, and by charge exchange with neutral grains: 

° So <=So S- So g° g° g- mTg_ at g_ ATg_ ^ go g_ go go ^ g_ go, gOgo v / 

a' a' a' 

We close this set of equations with constraints on the total number of grains in a given size bin, 

"g? + "g? + '^g^ = "g° ' (22) 
and the total number of an atomic species (neutral + positively-charged), 

"AO + — n-A , (23) 

and with charge neutrality: 

+ nj^+ - + X! ("-8+ ~ "g- ) ^ . (24) 

a 

This system of equations [( fTSl ) - (l24l) 1 is solved numerically via Gauss-Jordan elimination on the matrix equation derived by 
applying the Newton-Raphson iteration method. The rate coefficients in equations (fTSI l - (l24l) are given in AppendixlAl 

The mass of molecular ions is taken to be that of HCO+, — 29 m-p, while for the atomic ions an average value mp^+ — 
23.5 m-p, between the mass of Na {m^^+ = 23 Wp) and the mass of Mg (mMg+ = 24 rrip), is used. Since the ion masses are all 
comparable, the fact that different ionic species dominate in different density regimes does not affect the evolution of the cloud 
cores. The total number of atomic ions is fixed at 2.05 x 10^® ri„ (Morton 1974; Snow 1976). 

5. MAGNETIC FLUX LOSS AND ELECTRICAL RESISTIVITY 

The force equations [( [Tdl i - (fThl il and the induction equation (fTkl i are not written in the most convenient form for our purposes 
(see §|2]below). A useful simplification can be made, which amounts to a generalized version of Ohm's law; namely, we replace 
equations ( fTdb - ( fThl l with a modified form of equation ( fTkl ). This auspiciously eliminates 5 variables {v^, v-^, v^^, and 
Vgj^), but not without a cost. The ensuing algebra is messy, and much of it is deferred to Appendix lB.il Here, we outUne the 
simplification and highlight some results suitable for the present discussion. 



5.1. Resistivity of a Magnetic Gas 
The rate of change of magnetic flux across a given surface S, comoving with a fluid with velocity v, is given by 



dt Js[dt ^ ' 

Using Faraday's law (fTkl) . the integrand can be rewritten as 



■ dS . (25) 



dt 

and the current density can be calculated from 



= -c Vx [^E+-xBj -dS, (26) 

j = (T (i; + ^ X . (27) 

The quantity v is the velocity of the fluid, which for a weakly-ionized gas is essentially that of the neutrals Vn, and cr is the 
conductivity tensor. The presence of a magnetic field introduces an anisotropy in the equations, which is the reason for which the 
conductivity must be describ ed by a tenso r. If we take the 3 -direction to lie along the magnetic field, the conductivity tensor has 
the following representation (lParkslll99l]) : 

(0-_L -0-H \ 
au <J^ . (28) 
aj 

As B — > 0, the tensor must reduce to an isotropic form; i.e., cth — > and a± CT||. Because magnetic forces vanish along the 
magnetic field, cr| | must be independent of the magnetic field strength. 
Equation ( |27] i may be inverted to obtain the electric field, in which case a resistivity tensor rj is defined by 

E+'^xB^rij, (29) 

where, in the same representation as cr written above. 




\-m v± , (30) 



and 

V\\ = — , = ^5— — 2"' = — -■ (31a,b,c) 

The flux-freezing approximation corresponds to the limit rj ^ 
If we write the current density j in component form, it follows that we may write equation i 

£^ + ^ X B = Tallin + ri±j± + mj X b , (32) 

where j 1 1 and j j_ are the components of the current density parallel and perpendicular to the magnetic field, respectively, and 
b is a unit vector along the magnetic field. This relation between the electric field and the current density can be substituted in 
equation (|26] | to find that 

- ^ V X (?7|| j|| + r]^j^ + ijuj xb)-dS. (33) 



dt 



s 



This is the general form of the equation describing the loss of magnetic flux from a parcel of neutral gas, written entirely in terms 
of the components of the resistivity and current density tensors. 

For our model cloud, we have assumed axisymmetry and neglected rotation. In this case, the magnetic field is purely poloidal 
and the current density is purely toroidal by Ampere's law ([Ti]i. This geometry implies that the only nonvanishing component of 
the current density is the component perpendicular to the magnetic field, j = j±- The evolution of the poloidal magnetic flux in 
the neutrals' reference frame is then given by 

^ = -cJ^Vx{rj^j^)-dS. (34) 

The equivalent equation governing the evolution of the poloidal magnetic field is 

OB 

— = V X (-un X B) - cV X . (35) 

Equation (|34] | describes the evolution of magnetic flux in the neutrals' reference frame due to the motion of charges at right angles 
to the magnetic field and includes the effects of both ambipolar diffusion and Ohmic dissipation. The rate at which magnetic flux 
is lost equals the sum of the rates due to each process. Therefore, rj± can itself be written as the sum of two components, one 
related to ambipolar diffusion (subscript AD) and the other to Ohmic dissipation (subscript OD); 

VI- = V^AD + V-LOD ■ (36) 

The iss ue of how to separate the re s istivity r]± into its two componen t s is di scussed, for example, in iNakano & Umebavashil 
(ll986allbh . lGoldreich & Reisenegged (Il992l) . and lDesch & Mouschoviai(l2001h . We quote flie result here: 

V\\OD = V\\, V±OB=V\\, '7||AD=0, v±Au = V±-V\\- (37a,b,c,d) 
We now derive expressions for the resistivities from first principles. 
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5.2. Generalized Ohm 's Law 



We outline the derivation of a generalized Ohm's law, taking into account both elastic and inelastic collisions between neutrals, 
ions, electrons, and charged and neutral grains. We begin by writing down the force equation for the charged species s: 

= nsQs (e+ — xb)+ ^{Vn - Vs) + -^(t^go - Vs) . (38) 

The subscript s runs over all the charged species, taking on the values s = i, e, g+, and g_. Although we employ a grain size 
distribution, we consider only a single grain size in what follows for ease of presentation; a discussion of the consequences of a 
grain size distribution is deferred to Appendix IB. 21 The charge qg of species s carries an algebraic sign (e.g., it is negative for 
electrons). We write Tg inoi to represent the timescale for species s to be created by or take part in any inelastic collision. For 
example, 

+ + ^ + + (39) 

-)i,inol T'gog+jinel Pgn V'g+o.inol ''g+g_,inol ''g+go ,inel / . 

is the timescale for a neutral grain to participate in any inelastic reaction involving conversion between positive and neutral 
grains. The first two terms represent the production of positive grains due to charge exchange between neutral grains and ions 
and between neutral grains and positive grains, respectively; the next two terms represent the conversion of positive grains to 
neutral grains via neutralization with electrons and negative grains, respectively; and, the final term represents the conversion 
of positive grains to neutral grains via charge exchange. Since these are processes occurring in parallel, the reciprocals of their 
respective collision times are added to obtain the net collision time. Similarly, 

-1 

(40) 

."^gocanol 7"gog_,mcl Pgo V''"g_i,incl 7"g_g+,incl ''"g_ go,inol ^ 

is the timescale for a neutral grain to be involved in any inelastic reaction involving conversion between negative and neutral 
grains. The force equation for the neutral grains is 

= _ ^^J + ^ ^(^, _ , (41) 

where the index k runs over all the charged species. 

We eliminate the velocity Vg of species s in favor of a new velocity, Wg, which is the velocity of species s with respect to 
the neutral gas (Wg = Vg — v^)- In addition, we define as the electric field in the frame of reference of the neutral gas 
(En = E + Vn X B/c). Equations ( l38T l and dTTT l then become, respectively, 

= (^En + WgXb)- Wg + -^Wg„ , (42) 

1 + Qg \B J 1 + Qs 



= tOgo - 2^ Wk , (43) 

where we have introduced the cyclotron frequency of species 5, cj^ ~ qsB/msC, and have defined Qs and tq by 

Qg=P^^^ and 1 = ^+^^—. (44a,b) 

Equations ( |42] | and (|43] l form the set of equations to be solved.^ The species velocities (relative to the neutrals) Wg can be 
expressed in terms of E^ and then substituted in the definition of the current density 

3 = ^ngqgWg, (45) 

S 

where we have used charge neutrahty Ti^qs = 0). This expression can then be inverted to find in terms of j, which 
defines the resistivity tensor. The magnetic induction equation is then found by substitution into Faraday's law of induction: 

— - V X [vn xB)^ -cV X En . (46) 

Using this approach, we derive an induction equation generalized to include Ohmic dissipation, ambipolar diffusion, and the Hall 

effect for a six-fluid system including both elastic and inelastic collisions. 

The ensuing calculation is tedious, and we defer the details to Appendix lB.il Here we give only the final result: 

J = CT||£;n.|| +0-i-En,± - O-H-En X b, (47) 

where 



3 



The quantity (>s was written as in lTassis & MouschoviasI )2007al) . We have renamed it here to avoid confusion with the cylindrical radial coordinate 1 
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nsq^Tsn/iTis. The quantities Wg, fs, and are defined in AppendixiDl they 



The conductivity of species s is given by cr, 
represent the effects of inelastic collisions on the conductivity of the gas. In the absence of inelastic collisions, these formulae 
reduce to their standard form (e.g., see Parks 1991): 



s 

Equation (l47b may be inverted to give 



CTH 



GsUJsTau_ 
,2^2 



with the resistivities ?7||, 77 j^, and t^h given by equations dSTT i. 



(49) 



5.3. Attachment of Species to Magnetic Field Lines 

It is possible to write the velocity of each species, Vg, in terms of the velocity of the neutrals, v^, and the velocity of the field 
lines, , which is defined implicitly by 

(50) 



E + —XB = Q. 



The algebra and some intermediate results of interest are given in Appendix [Cj here, we quote the main result and explain it 
physically: 



Vs. I = V„ 



Vs X b ^ Vn X b 



1 



ViXb 



iV[-Vn) XbAs, 



1 



[Vf. 



I A, 



(51a) 
(51b) 



where the expressions for Qs and As are given in Appendix ICl The quantity 0^ is the attachment parameter (i.e., for Qs ^ 1, 
Vs K, Vi and species s is attached to the field hnes, wh ereas for Qs ^ 1, ~ I'n and species s is detached and comoves with 
the neutrals) — see, also. lCiolek & Mouschoviasl(ll993h . § 3.1.2. The function As quantifies the relation of one component of the 
species velocity to its mutually perpendicular component of the field line drift velocity, and essentially embodies Ampere's law. 
Under the assumptions of this paper, the midplane velocities of the charged species s, written in cylindrical coordinates (r, (f), z), 
are 

t)s,0(?-, Z = 0) = (Vn,r - ■Wf,r) As , Vs^r{r, Z ^ 0) ^ Vn,rT^ — + ,r ^ ! , ■ (52a,b) 



Os 



1 



Os 



1 



The first equation says that the charged species move in such a way as to cause differential motion between the field lines and the 
neutrals (Ampere's law). The second equation gives the radial velocity of any charged species in terms of the the velocities of the 
neutrals and of the field lines. These may be combined to yield 



Os V 



Os 



1 A, 



(53) 



In other words, the radial drift between species s and the neutrals is directly proportional to the contribution of species s to the 
azimuthal current. 

5.4. Grain Continuity Equation 
In the notation of Section ls!2l the grain continuity equation (fTbl ) may be written as 



^ + V.(pgt;„) = -V.(pg_ti;g 



where pg ~ pg_ + + pg^ is the total grain density. Eliminating Wg^ using equation ( |43] |. we find that 



dt 



V-ipgVn) = -V- 



Pg^Wg 



1 



To 



'g+n 



To 



'g-n 



(54) 



(55) 



We may use equation jCll l to eliminate the differential velocities of the charged grain species to find, after some manipulation, 

(56) 

where the components of the grain-continuity resistivity tensor, Jy^ont' ^i"^ defined as 



dp 

+ V • {pgVn) = -V • (77cont,||J|| + Vcont,±j ± + 7?cont,HJ X ^) 



f?cont,_L 



'?cont,|| 

E 



s=g+.: 



rus 



E 

s=g+,g- 



ms 



T0_ 



1 + —g, 

Ts 



T0_ 



1 + — ft 

Ts 



1 , 

1 H ft 



(57a) 
(57b) 
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T/cont,H 



E 



rus 



To 



These equations apply to all grain sizes separately. Under the assumptions in this work, j [ 
Equation ( l56b then becomes 



+ V • (pgVn) = -V • (?7cont,Hi X 



(57c) 

= and V • J = by axisymmetry. 

(58) 



Note that if 77co„t,H 
expected. 



0, a quantitative implementation of flux-freezing, the grain species are advected with the neutrals, as 



5.5. Joule Heating 

The rate Fdiff at which collisions dissipate kinetic energy as heat per unit volume (in the reference frame of the neutrals) may 
be calculated by taking the dot product of equation ( |42] | with Wg and using equation PTt : 



Tdiff — 




(59) 



where the summation indices s and k run, as usual, over all charged species (including charged grains of different sizes if a grain 
size distribution is considered). Using equation jClb to eliminate the velocities in favor of the current density, we find after much 
simplification 



Tdiff = V\\\3\ 



V±\3_ 




(60) 

0), this equation reduces to the usual 



In the limit where inelastic collisions are negligible relative to elastic collisions (i.e., Qs 
expression 

Tdiff '7||b'||P +??±b'±P ^ ??ODb'P + ?7ADb■±P• 
In the last step, we have used equations dSTT i to separate the contributions of Ohmic dissipation and ambipolar diffusion to the 
heating rate. Ohmic dissipation affects the total current density, whereas ambipolar diffusion affects only the perpendicular 
component of the current density. 

6. RADIATIVE TRANSFER 
6.1. The Flux-Limited Dijfusion Approximation 

Computing a formal solution of the full angle-frequency dependent non-LTE radiative transfer equation in a multidimensional 
numerical algorithm is a prohibitive task. Even if a rigorous yet tractable algorithm were developed to this end, the computational 
ex pense involv ed would prevent a solution in any reasonable amount of time. In fact, the sophisticated numerical code described 
in IStone et al.l (1992) designed to solve this problem with as few approximations as possible never saw public release. The 
ELD approximation is an attractive method for handling transport phenomena that is relatively easy to implement, robust, and 
inexpensive. It has the advantage over other diffusive approximations in that it preserves causality in regions where significant 
spatial variation can occur over distances smaller than a mean free path. For example, the Eddington approximation consists of 
assuming the radiati on field is everywhere iso tropic, an assumption that is violated in the optically-thin limit where the radiation 
becomes streaming f Mihalas & Mihalaill984h . 

The fundamental assumption of ELD is that the specific intensity is a slowly varying function of space and time. This is certainly 
valid in the diffusion and streaming limits (at least in one d imension); one hopes tha t it is approximately true in intermediate 
situations (and in multi-dimensions). Given this assumption. iLevermore & Pomraningl(ll981l) showed that the radiation flux can 
be expressed in the form of Eick's law of diffusion, 

^ = -I?FLD , (61) 

where the diffusion coefficient I?fld can be written as 

cAfld 



(62) 



The dimensionless function Afld ~ Afld(^^) is called the flux limiter. Similarly, in ELD theory the radiation pressure tensor 
can be expressed in terms of the radiation energy density via 

P = f£ , (63) 
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Fig. 2. — (a) Planck mean absorption coefficients and (b) Rosseland mean extinction coefficients (in cm-^ per gram of dust) for grain sizes a = 0.0256 /xm 
(solid line), 0.0543 /im (dotted line), 0.1190 fim (dashed line), 0.2600 fim (dash-dot line), and 0.5680 /im (dash-dot-dot-dot line). 

where the components of the Eddington tensor f are given by 



f =^(l-/)l + ^(3/-l)nn, 



(64) 



where n = V£/|V£ is the normahzed gradient of £ and the dimensionless function f — f {£) is called the Eddington factor 
dTurner & Stonell2001h . The flux limiter Afld and Eddington factor / are related through implicit constraints between the 
moments ^ and P, so that 

/ = Afld + ^Ild"^^ > (65) 

where TZ is the dimensionless quantity TZ = \'V£\/xj^£- We have chosen the flux limiter derived by Levermore 8c Pomraning 
(1981, eq. 28), which is given by 

2 + n 

- 6T3^W ■ ^^^^ 

Its use in hydrodynamic simu lation s of star formation has b een documented, for example, in iBodenheimer et all (119901) . 



Its use in nygroaynamic simu lation s or star lormation nas p ei 
IBodenheimer et alJ(ll993Lll995h . and lWhitehouse & Batd (l2006h . 



6.2. Dust Opacities 

For temperature less than ~ 1500 K, the contribution of dust to the total opacity dominates that from all other sources. We take 
ke = Kp and x:f — XR (see Mihalas & Mihalas 1984, § 82), where Kp and xr have been obtained from private communication 
with Dmitry Semenov and Thomas Henning . The major dust co nstituents are "iron-poor" silicates, troilite, organics, and water 
Their relative mass fractions are taken from Pollack et al.l (|1994|) . These opacities (in cm^ per gram of dust) are shown in Figure 
|2]for the five different grain size bins taken to represent an MRN distribution (see § 4.2 above). The major changes in the dust 
opacities are: for temperatures T < 120 K, all dust material are present; at T ~ 120 K, water ice evaporates; at T = 275 K, 
volatile organics evaporate; at T = 450 K, refractory organics evaporate; at T = 680 K, troilite (FeS) evaporates. 

7. MODIFIED ZEUS-MP CODE 

In order to solve for the evolution of the many complex, nonlinear systems of equations presented in this paper, numerical 
techniques are necessary. Rather than extend the sophisticated fully implicit, nonorthogonal adaptive mesh, two-fluid MHD code 
of Fielder & Mouschovias (1992) to include radiative transfer, we have opted instead to modify the publicly-available Zeus- 
MP RMHD code dHaves et al.ll2006l) . We have altered the algorithms governing the evolution of the magnetic field in order to 
account for ambipolar diffusion and Ohmic dissipation. In addition, we have added routines to evolve the total grain density and 
to compute the species abundances from the equilibrium chemical model detailed in § 14.31 and have made changes to Zeus-MP's 
adaptive mesh module in order to track the collapsing core. New modules were also written to improve both the efficiency of 
Zeus-MP's implicit radiative transfer solver and the manner in which the gravitational potential is calculated. A brief description 
of these modifications follows. 

The evolution of the magnetic field is governed by equation ( [35] l. This equation does not assume flux-freezing in any species and 
accounts for both ambipolar diffusion and Ohmic dissipation. The first term on the right-hand side, V X {v^ X B), represents 
the advection of the magnetic field by the neutrals. As long as the nonideal MHD Courant condition 



At < At 



difi' 



4tt {Axf 

C^TJ± 2 



(67) 



is sat isfied, then the method of characteristics used to update the magnetic field due to this term remains valid jMac Low et alj 
Il995h and we may use it without modification. Physically, this is because, for a sufficiently short timestep, ambipolar diffusion 



14 



does not have time to alter the characteristics of Alfven waves. The second term on the right-hand side, V X [(c^77j^/47r) V X B], 
represents the diffusion of the magnetic field lines and is applied using the method of constrained transport. This is done in the 
source step part of the code. We employ a similar approach to the grain continuity equation ( ISSl l. The right-hand side is treated 
as a source term in the source step part of the code. Then, the grain mass density is advected during the transport step using the 
multi-species advection module already present in Zeus-MP. Joule heating (eq. Il60l ) is also applied to the internal energy during 
the source step. 

We use an adaptive grid to track the evolution of the contracting core. The grid, which must resolve the core, has its innermost 
zone constrained so as to always have a width in the range AT,cr/5 — At,ci/10, where At.ci = lACrg is the critical thermal 
lengthscale (iMouschoviasll 1 99ll) . C = ik-QT /my^Y^'^ is the isothermal sound speed of the gas, and Tff = (37r/32Gpn c)^^^ is 
the spherical free-fall time. (The quantity ^ is the neutral mass density at the cloud center). The critical thermal lengthscale 
is the smallest scale on which there can be spatial structure in the density without thermal-pressure forces smoothing it out. The 
number of cells is fixed, and their positions are spaced logarithmically, so that the spacing between the ith and ith + 1 cells is a 
number greater than the spacing between the ith — 1 and ith cells. 

As discussed in § 3.8.1 of Haves et al. (2006), Zeus-MP uses the diagonal preconditioned conjugate gradient method to solve 
the sparse matrix equation that results from spatially discretizing equations ( llml i and ( fTnl i. Diagonal preconditioning is an at- 
tractive technique due to its simple calculation, the fact that it poses no barrier to parallel implementation, and its fairly common 
occurrence in linear systems. However, it is only efficient for matrices in which the main diagonal elements are much greater in 
magnitude than the off-diagonal elements (a condition referred to as "diagonal dominance"). Unfortunately, this is generally not 
the case for the problem studied here. We have therefore replaced the diagonal preconditioner with an incomplete Cholesky de- 
composition preconditioner, similar to what was provided in the public-release version of Zeus-2D. The savings in computational 
cost has been enormous. 

Zeus-MP computes the gravitational potential through a two-step process: first, the gravitational potential ijj is found on the 
computational boundaries; then ijj is found in the interior by iteratively solving the Poisson equation for ijj using a sparse matrix 
solver that relies on the preconditioned conjugate gradient method. For each boundary of the domain, there are two possible 
boundary types: (1) Neumann, in which the slope of the gravitational potential is set to zero; and, (2) Dirichlet, in which the value 
of if} in the ghost zones is specified. (There is actually a third possible boundary type — periodic boundary conditions — however, 
this boundary condition is not used here.) Neumann boundary conditions are used at symmetry boundaries (axis r = and 
equatorial plane z — 0), while Dirichlet conditions are applied at the outer boundaries, far from most of the mass distribution. In 
the public -release version of Zeus-MP, Dirichlet boundary conditions are implemented by computing -tp on the domain boundaries 
using a multipole expansion formula, which we give here in spherical coordinates ( ^, 9, (j)) for an axisymmetric mass distribution: 




where P( is the Legendre polynomial of degree i. Note that Zeus-MP uses only the monopole and quadrupole moments, in 
contrast to the earlier Zeus-2D code (Stone & Norman 1992), which used arbitrarily high £ moments until a desired convergence 
was achieved. The term in brackets is denoted by qi, and is known as the multipole moment of order £ of the density distribution. 
In most situations, a dozen or so multipole moments are sufficient for convergence. They can be calculated once and then 
used to find the potential at many boundary points. However, we have found this subroutine inadequate for our purposes, as 
it fails to converge in situations when the distance ^^to the point at which one wishes to calculate the potential is greater than 
the distance to any mass element. In the axisymmetric geometry used in this work, this situation is inevitable: mass elements 
near (r, z) = (0, Z) or (r, z) = (R, 0) are closer to the origin than are mass elements near (r, z) = (R, Z). We therefore have 
followed iDesch & Mouschoviasl(l2001h in using the more general multipole expansion (lJacksonlll999h : 

oo 

where («;<) is the greater (lesser) of ^' and ^, and £ may take arbitrarily large integral values until a desired convergence is 
achieved. An unfortunate consequence of this more general expansion is that it is not possible to perform one integration over all 
space and use the result of that integration (the multipole moment qi) to find the potential at all boundary points. Instead, a new 
integration over space must be performed for each boundary point, since the location of that boundary point will determine how 
the integral in equation f69[ is separated into two integrals: one integral will have ^' in the numerator of the integrand, and in the 
other integral ^' will be in the denominator. This situation is further complicated by the use of parallelization, since comparisons 
of and and subsequent integrations must take place across multiple processors. 

8. SUMMARY 

We have formulated the problem of the formation and evolution of fragments (or cores) in magnetically-supported, self- 
gravitating molecular clouds in two spatial dimensions. The evolution is governed by the six-fluid RMHD equations. The 
magnetic flux is not assumed to be frozen in any of the charged species. Its evolution is determined by a newly-derived gener- 
alized Ohm's law, which accounts for the contributions of both elastic and inelastic collisions to ambipolar diffusion and Ohmic 
dissipation. The species abundances (electrons, atomic and molecular ions, positively-charged grains, negatively-charged grains, 
and neutral grains) are calculated using an extensive equilibrium chemical network. Both MRN and uniform grain size distribu- 
tions are considered. The thermal evolution of the protostellar core and its effect on the dynamics are followed by employing the 



p{i.',9')Pe{cose')- 



i+i 



Peicosd) , 



(69) 
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grey FLD approximation. Realistic temperature-dependent grain opacities are used that account for a variety of grain compo- 
sitions. We have augmented the publicly-available Zeus-MP code to take into consideration all these effects and have modified 
several of its algorithms to increase its accuracy and efficiency. 

We summarize here for convenience the simplified evolutionary equations discussed above and used in our modified version of 
the Zeus-MP code: 

■ V.(p„t;„) =0, (70a) 



dt 

dt 



V • (pgi;„) = -V • (?7cont,Hi X . (70b) 



' V • {p^v^v^) = -VP„ - pV^ + -^(V X i?) X S - AfldV^ , (70c) 



— = V X Un X B - — ^ V xB] , (70d) 



V^tl> = inGpn , (70e) 



dB ^ ( ^ c^Tj 
— = Vx [vnXB- — 

— ^ + V • (Unt^n) = -^^n V • - AnKpB + CKpS , (70f) 

ot 

— + V • (£Vn) = V • ( ^-^Vf I - VVn : P + ilTKpB - CKpS . (70g) 
ot V XR / 

These equations are considered together with the relations Pn = (7— l)unandP — f£. Results will be presented in a forthcoming 
paper (Kunz & Mouschovias 2009, in preparation). 
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John Hayes for his assistance with the Zeus-MP code; and, Dmitry Semenov and Thomas Henning for generously providing the 
dust opacities. TM acknowledges partial support from the National Science Foundation under grant NSF AST-07-09206. 



APPENDIX 
A. RATE COEFFICIENTS 

For radiative recombination of atomic ions and electrons, a^- = 2.8 x 10~^^ (300 K/T)°-^^ cm'^ s~^; for the dissociative 
recombination of electrons and HCO+ ions, adr = 2.0 x 10"'' (300 K/T)" "^^ cm^ s"^ (lUmebayashi & Nakariol 1990). The rate 
coefficient adopted for charge exchange reactions between atomic and molecul ar ions is (3 = 2.5 x 10 ~^ cm'^ s~^ (Watson 1976). 

The rate coefficient s involving gas-phase species and grains are taken from lSpitzed (Il941lll948l) . with refinements made by 
iDraine & SutinI (Il987h to account for the polarization of grains: 
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The sticking probab ilities of electrons or ions onto grains, denoted Vc and Vu are assigned the values 0.6 and 1.0, respectively 
(lUmebay ashill 1 9 8 31) . Other quantities in these equations have their usual meanings. 
The rate coefficients for charge transfer between charged grains are given by 



Tra, 



SkBT \ 
Trmrcd / 



1/2 



r skBT 



1/2 



1 



2 + {a, 
.2 \ 1/2' 



2as 



Vaa' , 



1/2' 



(A5) 



(A6) 



where the reduced mass of the two grains (labeled a and a') is defined by mrcd — mama' / {ma + ma'), and Osum = + Oa' is 
the sum of the radii of two grains a and a'. The probability of two oppositely charged grains neutralizing each other upon contact 
is assumed to be unity. The probability of charge being transferred to a neutral grain gp from a charged grain g^ is assumed to 



be proportional to the surface areas of the grains, so that Vaa' — '^aa'/i'^a 



In other words, of all the collisions between a 



neutral grain, a, and a charged grain, a', only a fraction Va lead to charge exchange. The complementary probability Va' leaves 
the charges unchanged. 
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B. GENERALIZED OHM' SLAW 
B.L Derivation 

We consider equations (l42T i and ( |43T l. repeated here for convenience, which are to be solved for the species drift velocities Wg 
relative to the neutrals: 







— En+Ws Xb] -W 



Qs 



l + Qs 



= -f^go - XI 



To 



-Wk 



(Bl) 
(B2) 



We define, for brevity and clarity of presentation, the following quantities 

To LOsTs-a 



Tsancl 1 + gls 



1-E 



To gfc 
Tfc,inel 1 + Qk 



To ^sTsn 
Ts.inol (1 + eis)^ 



To gfc 
Tfcanol 1 + Qk 



To UJsTs-a Qs 

Tsancl 1 + gs 1 + 
To Qk 

Tfc.incl 1 + gfc 



1-E 



*2 = E*2,^ 



^3 = E *3,fc ■ 



(B3a) 



(B3b) 



(B3c) 



We recall that the index k runs over all the charged species independently of the index s, which denotes the charged species in 
question. Note that the denominator in the above expressions may be written with the help of equation (l44b) as 



To 



1 



'gon 



Tfe,incl 1 + gfc ' 



which shows its positive definite nature. 

We first multiply equation (IB 11 1 by tq/ts inci, sum over s, and use equation ( IB2l i to find that 



.„ = ^i^E,, + ^^i,kWk X b, 



(B4) 



(B5) 



where we have switched the summation index to k to avoid confusion with the species in question, s. Next we take the cross 
product of equations ( IBU and ( IB5I ) and the unit vector b: 



w, xb = 



1 + gs vs 



1 + g. 



■w^„ X b, 



w„ 



go X b = fl-^^n X 6 - X *l,fcU^fc, 

k 

Equation ( |B7| i is now substituted in equation (|B6j to obtain 



tfJs X 5 = 



l + gs ^ 



1 + gs 1 + gs 

Inserting this expression in equation (IB5l l. we find that 

Ifgo = *1-|-E„ + (*2 + *3*l)-|£^n X 6 - y ('J'2,fc + 4'3*l,fc)lf'fc,± 



(B6) 
(B7) 

(B8) 
(B9) 



which is now ready to be inserted, along with equation ( IB8b . in equation (IBU : 
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The symbol Ssk is the Kronecker delta. This is our first main result: it gives the velocity of each charged species in terms of the 
electric field in the frame of the neutrals. Another way of interpreting this equation is obtained by defining the velocity of the 
magnetic field lines with respect to the lab frame: 

V{=-^Exb. (Bll) 

Then equation ( IB 10b provides the velocities of all the charged species in terms of the neutral velocity and the field-line velocity. 
We made use of this concept earlier in Section |53] 

Equation ( IB 10b can be separated into components parallel and perpendicular to the magnetic field. The parallel component of 
the current density is easily obtained: 

j|| = ^n^QsW^^ii (B12a) 

S 

= !:•.* (ff^ + Yf^*.)|E,„N (B.2b) 

= ^a,(l-<;,)£;n^ll (B12c) 

s 

= a||£;„.l|, (B12d) 

where we have introduced the conductivity of species s, as — nsq^Tsn/nis, and c^s, given in Appendix iDl is the factor by which 
the conductivity of species s is altered because of inelastic collisions. In the last step above, we have introduced the parallel 
conductivity, (t^\, which is defined in situ. Note that <;s > for all s. In other words, by interfering with the rate at which 
the charge carriers flow along the magnetic field, inelastic collisions are responsible for decreasing (increasing) the parallel 
conductivity (resistivity) of the gas. 

Finding the perpendicular components of the current density is not as straightforward and amounts to solving a matrix equation. 
We first define the 4x1 column vectors and C^, whose entries are given by 

(C-^)^ =a;,T,„(l-<j,) and (C")^ = - • (|B]13a,b) 

We also define the 4x4 matrix of coefficients A whose entries are given by 

(A)^^ = [1 + o-^rfji _ ^^)] + u;lTl^sk{l - Ssk) • (B14) 

The expressions for <js, zJs, 'Ps, and dsk are given below in Appendix ID] Then the perpendicular component of equation dBlOb 
takes on the form 

C^^E,,,^-C^^E,,xb = AW^, (B15) 

13 13 

where W± is the 4x1 column vector of unknown velocities of charge species relative to neutrals, [w^, w^, Wg_ , lOg^]''". 
We use Cramer's method to solve the matrix equation ( |B15b . We define 

D = det[A]. (B16) 

In addition, we use the notation Djr to represent the determinant of A with the sth column of A having been replaced by C^. 
Similarly, is the determinant of A with the sth column having been replaced by C^. Then, the solution of the system ( |B15b 
is 

c c 

"^'■^ " "d" fi-^"'^ " "d" B ^" ^ • ^^^^^ 
Once the determinants have been computed, the current density perpendicular to the magnetic field may be obtained: 

j± = ^ nsqsWs,± (B18a) 

s 

^ ^ " En _L - " EnXb (B18b) 

D B ' D B 

- E TT^w^"^(^)^-- + E ;r?"i!r"^V M^„ X b (B18C) 

= a±En,± - aiiEn X b . (BlSd) 

In the last step, we have defined the perpendicular conductivity a± and the Hall conductivity (Th, which include the effects of 
inelastic collisions. 
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B.2. Modification due to a Grain Size Distribution 

When considering a grain size distribution, rather than single-size grains, two changes must be made to the above derivation. 
First, the inelastic collision timescales given by equations (|39] l and ( l40l i must be modified as follows: 
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(B19a) 



(B19b) 



The summations over a' (the grain size label) indicate that a grain of size a may give or receive charges not only from other 
grains of its own size, but also from all other different-size grains. Second, the summation index s in equations (IB12l l and ( |B18l l 
should range over all charged species, including all sizes of charged grains. 

C. DERIVATION OF SPECIES VELOCITIES 

En route to the derivation of a generaUzed Ohm's law, the differential velocity of every species can be obtained in terms of the 
current density: 

= <^\\,sV\\j\\ + <^±.s{v±3± + mo X 6) - crH.s(T?±j X 6 - mJA) 

Using equation jBlll l. is is straightforward to show that 



C77H . 

W{ = Vf -Vn = -^3 Xb —^3i_ ■ 
We may then write the components of the current density in terms of the differential velocity of the field lines as 
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(TnWf ± + a±Wf X b . 



Defining the indirect coupling coefficient implicitly by 
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A,, 



B 



B 



cusqs 



equation (IClb may now be written in component form as 
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or, more explicitly. 
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These equations were discussed in Section 



D. DEFINITIONS 



In the main text, as well as in the preceding appendices, we had delayed giving explicit definitions of <js, tUj, (ps, T"s, and iDsk 
for all (s, k) = e, i, g_, and g+ due to their complexity and length. Here we give explicit expressions for these quantities for all the 
charged species. Before we proceed, however, a few simplifications are in order. Since both (rc^nei/Ten) and (ri^nei/Tin) ^ 1 
for the density regi me of interest in this paper, we may neglect the influence of inelastic collisions on the electron and ion fluids. 
Using the results of lTassis & Mouschoviasl(l2007bh . we also may assume that the velocity difference between a given grain and a 
neutral particle is less than the sound speed of the gas. These are both excellent assumptions, and lead to a much more compact 
form of the following definitions than would otherwise be possible. 

The variable firs t app eared in the definition of the parallel conductivity (IB 12b and again later in the definition of the perpen- 
dicular conductivity (IB18b . For electrons and ions, <;o = Q = 0, because of the negligible influence of inelastic collisions on the 
electron and ion fluids relative to that of elastic collisions. The expressions for the negative and positive grains are given by 
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and are clearly positive. As mentioned in Section 15.21 by interfering with the rate at which the charge carriers flow along the 
magnetic field, inelastic collisions are responsible for decreasing (increasing) the parallel conductivity (resistivity) of the gas. 

The derivation of the perpendicular conductivity involved many more definitions, all of which are given below. For the same 
reason stated above for which — Q — 0, the expressions for Ws and ips vanish when s = e or i. The quantity Tg is equal to 
unity for these species. The nontrivial rUg, fs, and Tg for s = g_, g_|_ are given by 
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In equation (IB 141) . we had introduced -dgk as a measure of the inelastic collisional coupling between different pairs of charged 
species, s ^ k. This variable was also used in the definition of Tg above. Since the effect of inelastic collisions on the electron 
and ion fluids is negligible, dgk vanishes for (s, k) = e or i. The only remaining nonzero values of dgk involve the charged grain 
species, and are given by 
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